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Abstract 

We propose a stochastic individual-based model of the progression of atrial fibrillation (AF). The model 
operates at patient level over a lifetime and is based on elements of the physiology and biophysics of AF, making 
contact with existing mechanistic models. The outputs of the model are times when the patient is in normal 
rhythm and AF, and we carry out a population-level analysis of the statistics of disease progression. While 
the model is stylised at present and not directly predictive, future improvements are proposed to tighten the 
gap between existing mechanistic models of AF, and epidemiological data, with a view towards model-based 
personalised medicine. 
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1 Introduction 

Atrial Fibrillation (AF) is an electrical heart disease 
characterised by irregular heart rhythm and absence of 
organised atrial activity, with a prevalence of 1 — 2% 
in the UK [1]. Incidence and progression of AF in¬ 
creases with age [1], and is correlated with hyper¬ 
tension, diabetes and other cardiovascular conditions 
[2, 3, 4]. Furthermore, AF significantly increases risk 
of stroke [5, 6, 7]; thus, the importance of understand¬ 
ing this disease and of predicting its progression is 
undisputed. Aside from improved monitoring of pa¬ 
tients and clinical intervention, computational mod¬ 
elling is becoming increasingly valuable in the medi¬ 
cal sciences, in particular with a view towards person¬ 
alised medicine [8]. A natural question is then to as¬ 
sess what is the most appropriate time and space scale 
over which to construct models of AF. Biophysically 
detailed models of cardiac cells and tissues can provide 
insights into the mechanisms that initiate and sustain 
AF [9], whereas population-level models, parametrised 
by available epidemiological data enable treatment 
pathways to be evaluated [10]. The most appropriate 
model will of course depend on the specific problem in 
hand. 

Biophysically detailed models of cardiac cells and 
tissue are based on the physical laws governing the 
electro-mechanics of the heart. Simulating these mod¬ 
els however is computationally expensive, and so only 


time scales of a few heartbeats in a single patient are 
accessible [8]. This constraint prevents a bottom-up 
simulation of progression of the disease in a single pa¬ 
tient at a time scale longer than a few minutes, and 
so these models are not suitable for complete study of 
the self-reinforcing effect of AF (‘AF begets AF’ [11]) 
which occurs over days, weeks, and years. Construct¬ 
ing predictive models for the life-time progression of 
AF in individual patients based on the underlying bio¬ 
physics is therefore difficult given current computa¬ 
tional resources. 

Population-level models can capture long time 
scales, but they require epidemiological data as an 
input to determine model parameters. Clinical data 
are available at medium time scales of several days, 
but these data are sparely sampled, and mostly col¬ 
lected from patients with implanted monitoring de¬ 
vices. These patients are more likely to have other 
cardiac-related co-morbidities, and so the sample is 
biased. In addition, the storage capacities of devices 
limit the temporal resolution of the data. To our 
knowledge, on a population-level, a quantitative model 
describing the progression of the disease at a long time 
scale is not yet available, and so the gap between acute 
mechanistic approaches and modelling on time scales 
of the human lifetime remains. 

The aim of the present study is to make progress 
towards closing this gap and linking the small scale 
biophysical models with epidemiological models at a 
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population level. Ultimately, one may hope that one 
day fully developed multi-scale models will be avail¬ 
able, capable of predicting the progression of a given 
patient through the stages of AF over a time scale of 
years and starting from personalised parameters (e.g 
extent of atrial fibrosis obtained from cardiac imag¬ 
ing). Such models, if they can be constructed, would 
provide a key step towards computer-based predictive 
personalised medicine, and the ability to more accu¬ 
rately design individual treatment plans and to predict 
the cost-effectiveness of new interventions population¬ 
wide. 

To work towards this long-term goal we propose 
a coarse-grained patient-level stochastic model, based 
on some elements of the current understanding of the 
mechanisms that spontaneously initiate, sustain, and 
terminate AF. The finest scale described by the model 
is the status of a patient with respect to the disease, 
i.e. a patient can at any one time be in an episode of 
AF, or in sinus rhythm (SR). The switching between 
these states is governed by two rates, describing initi¬ 
ation and termination of episodes respectively. Both 
of these rates are patient-specific in the model, and 
reflect genetic disposition, age and individual patient 
history. Most crucially, the model describes the ef¬ 
fects of episode-based remodelling of the heart, mak¬ 
ing future episodes more likely. To capture this AF- 
begets-AF effect the model retains a transient mem¬ 
ory of earlier episodes, and these affect the initiation 
and termination rates of later episodes. Many of these 
processes and the parameters describing them are— 
in principle—linked to the underlying biophysics. We 
stress clearly that establishing this link quantitatively 
is of course a major enterprise, and not within the 
scope of the present work. The spirit of our paper is 
very much to make a first step in this direction. 

The remainder of the paper is structured as follows: 
In the next section we describe some of the main fea¬ 
tures of AF, both from the biophysical microscopic 
perspective and from the long-term epidemiological 
view. Our account is not comprehensive; we focus 
on features relevant to our modelling approach. In 
Section 3 we then present the actual mechanics of the 
model with a detailed quantitative description avail¬ 
able in the Appendix for clarity. The outcome of our 
simulations is then described in Section 4, before we 
summarise our model and discuss its potentials and 
limitations. 

2 Main features of AF relevant to the 
model 

2.1 Biophysical processes 

We briefly summarise the biophysical phenomena con¬ 
tributing to the key processes described by the model: 
initiation and termination of episodes of AF, ignoring 
external interventions such as pacemakers, cardiover¬ 
sion, or effects of medication. The likelihood of initia¬ 
tion and termination of AF in the model is dependent 
upon remodelling of structural and and functional el¬ 
ements of the heart, which occur over multiple time 
scales [13, 14, 15]. 


2.1.1 Initiation of episodes of AF 

The transition from sinus rhythm (SR) to AF is char¬ 
acterised by a combination of a triggering event and 
an existing vulnerability of the atria to such trig¬ 
gers [16]. This combination results in self-sustaining 
rotating waves of electrical activation, which act to 
suppress normal beats arising from the heart’s natu¬ 
ral pacemaker. Examples of triggering events are a 
rapid increase of the atrial rate, and ectopic electrical 
beats within the atria [17, 18]. Increased vulnerabil¬ 
ity or triggers can be due to structural abnormalities 
in the atria such as fibrosis, hypertrophy or gap junc¬ 
tional lateralisation [19, 20, 21], and specific electrical 
changes, which transiently modify the electrical be¬ 
haviour of cardiac cells [3, 11, 16, 22, 23, 24]. 

We adopted a coarse-grained picture of these de¬ 
tailed electrophysiological processes underpinning the 
initiation of episodes of AF. In our model the physio¬ 
logical state of the atria of a patient is described by a 
binary state: the patient can either be in AF or SR. In 
our approach, we separate triggering events from the 
vulnerability of the atria. Trigger events are seen as 
external to the model, so we do not include changes in 
patient-specific behaviour or physiology affecting the 
frequency of triggers. For simplicity, we assume that 
triggering events occur at a constant rate, throughout 
the life time of a patient. The time between two trig¬ 
ger events is random and exponentially distributed. In 
our model the vulnerability then describes the prob¬ 
ability that a given trigger event leads to an episode 
of AF: the higher the vulnerability of the atria, the 
greater the probability that a triggering event induces 
an AF episode. 

We consider three factors contributing to a patient’s 
vulnerability: (i) a genetic predisposition, (ii) an age 
and co-morbidity related factor, and (iii) previous his¬ 
tory of the patient. These factors were succinctly sum¬ 
marised by Heijman et al. [12]. Based on this existing 
work we illustrate the different components contribut¬ 
ing to vulnerability in Fig. 1. Genetic pre-disposition 
is constant and does not depend on age. The contribu¬ 
tion related to age and co-morbidities increases with 
time; we assume implicit correlation of age to other co¬ 
morbidities such as pulmonary hypertension, hyper¬ 
trophy etc, and do not therefore consider independent 
effects of specific co-morbidities in the present study. 
The third element—previous history of the patient— 
models the memory effects of AF. The vulnerability of 
a patient at any one time depends on the time spent 
in AF previously and on the number of episodes the 
individual has experienced so far, and when these oc¬ 
curred. 

Factors (ii) and (iii) capture the effects of physiolog¬ 
ical remodelling: Structural and electrical properties 
of the atria are modulated by risk factors including 
ageing and previous history of AF [3, 16]. All factors 
are patient-specific: different patients in the simula¬ 
tion can, in principle, have distinct parameter settings. 

2.1.2 Recovery 

The biophysical mechanisms leading to spontaneous 
termination of an episode of AF are much less studied 
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Figure 1: A framework of AF initiation, maintenance and progression, modified from [12]. AF onset is dependent 
on vulnerability and triggers. Vulnerability is dependent on a constant genetic predisposition, time-varying age/co¬ 
morbidity-related remodelling and AF-related remodelling dependent on AF history. An increase in vulnerability 
may allow initiation of AF episodes upon a trigger event occurring; over time, some patients progress on to 
paroxysmal, persistent and permanent AF as their AF vulnerability increases. Note the timescale for AF triggers 
and episodes is distinct from the lower axis and is expanded for visualisation. 


and understood than those leading to the initiation of 
an episode. Our formulation is based on the hypothe¬ 
sis that the factors contributing to the rate with which 
episodes terminate are similar to those contributing to 
activation. Any factor leading to an increase in the ini¬ 
tiation rate is taken to cause a reduction in the termi¬ 
nation rate. Specifically, we assume one contributing 
factor reflecting age and co-morbidities. This term in 
the recovery rate decays (exponentially) as the patient 
ages, and does not depend on previous history of AF 
progression. A second component depends on previ¬ 
ous history: the propensity to recover from an episode 
of AF is reduced by time spent previously in the AF 
state. 


2.1.3 Acute clustering 

Recurrence of an AF episode shortly after an exist¬ 
ing episode has been observed in clinical studies [25]. 
Kaemmerer et al [26] suggest that the likelihood of re¬ 
currence is highest immediately after termination and 
decreases over time, and that a long episode is likely to 
be followed by multiple shorter episodes immediately 
after. In order to capture this effect we assume that ac¬ 
tivation rate is increased immediately after an episode 
and exponentially decays to baseline, which is inde¬ 
pendent of AF history. We assume the recovery rate 
follows a similar trend; however we already assumed 
that the history component of AF recovery rate tran¬ 
siently decreases as time in an episode increases. To 
capture the effect that long episodes are immediately 
followed by short episodes, we thus assume the recov¬ 
ery rate ‘overshoots’ above baseline immediately after 
an episode terminates, then decays back to baseline. 


2.2 Progressive stages of AF 

Three stages have been broadly defined clinically in the 
progression of AF: (i) paroxysmal AF , (ii) persistent 
AF , and (iii) permanent AF. This is briefly illustrated 
in the schematic of Fig. 1). 

In paroxysmal AF, a fibrillation episode has a typ¬ 
ical duration of minutes or hours, lasting up to 7 
days, with episodes spontaneously terminating. Inter¬ 
episode times can vary greatly from days to months 
[27], and inter-episode times will often reduce as the 
condition progresses, due to physiological remodelling 
[11, 24, 28]. The next stage, persistent AF , is indi¬ 
cated when the patient has AF episodes lasting 7 days 
or more, which do not self terminate and require med¬ 
ical intervention such as cardioversion. A patient is in 
permanent AF when the episodes have continued for 
more than 1 year; typically interventions to terminate 
the fibrillation will fail. 

We stress that although these stages are used to se¬ 
lect appropriate treatment strategies, it is difficult to 
diagnose the stage any one AF patient is in clinically. 
Any such assessment is therefore approximate in prac¬ 
tice. 

2.3 Model Outline 

Our model describes a patient as either in AF or in 
sinus rhythm at any one time. A patient trajectory is 
therefore a series of points in time at which episodes 
begin and end. The starting point of the simulation is 
the birth of the patient, and all episodes are ‘detected’ 
in the simulation. A real patient would not necessar¬ 
ily be diagnosed to have paroxysmal AF upon initial 
occurrence of an AF episode; thus we use the follow¬ 
ing classification in our simulations, (i) A patient is 
considered to progress to paroxysmal AF immediately 
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(a) Main elements governing the switching process be- (b) Time courses of vulnerability and recovery 
tween sinus rhythm and AF. rate during AF episodes and in sinus rhythm 


Figure 2: Illustration of the main mechanics of the model. Vulnerability of AF due to ageing, V age , increases 
according to a sigmoid function, and age-based recovery rate Rq decreases exponentially. A healthy patient 
triggers an AF episodes with rate V(t)/r and terminates an AF episode with rate R(t). In AF, episode-based 
remodelling V ep i transiently increases vulnerability at a rate a and decreases the total recovery rate R at rate fi. 
When an episode terminates, the recovery rate ‘boosts’ above baseline by a factor of B. Back in sinus rhythm, 
Wpi decays to zero with rate f3 and R ep i decays with rate v back to the age-based baseline Rq. 


following an episode of at least 14.4 minutes (1% of 
a day), (ii) They are said to be in persistent AF fol¬ 
lowing an episode lasting 7 days, and (iii) progress to 
permanent AF from the point in time at which they 
enter an AF episode which continues until the end of 
the simulation at age 100 years, without return to si¬ 
nus rhythm. 

3 Description of the model 
3.1 General setup 

The model describes the trajectory of individual pa¬ 
tients throughout their life time. Patients can be sim¬ 
ulated individually, and they do not interact. Time in 
the model is continuous, and at each point in time t a 
patient is characterised by three attributes, their AF 
status S(t), their vulnerability, V(t), and their recov¬ 
ery rate R(t). 

The AF status indicates if the patient is in AF at 
time £; if so we write S(t) = 1, otherwise S(t) = 0. 
As far as the occurrence of AF episodes is concerned a 
patient trajectory is fully characterised by the values 
S(t) for all times t. Simulations start at age of zero 
years, this is labelled t = 0. They end at age 100 years 
(t = 100). 

The behaviour of S(t) is determined by the activa¬ 
tion and termination rates of AF episodes. If a pa¬ 
tient is in sinus rhythm at time t they can begin a 
new AF episode, this occurs with rate V(t)/r in our 
simulation. The prefactor r _1 indicates the frequency 
of trigger events. As described in earlier sections the 
initiation of episodes of AF occurs based on trigger¬ 
ing events combined with existing vulnerability of the 
patient. The factor r _1 reflects the trigger, and V(t) 
stands for the vulnerability. 

If a patient is in an episode of AF at time t ( S(t ) = 
1) they may return to sinus rhythm. This occurs with 
rate R(t). 

Mathematically speaking the model describes a bi¬ 


nary process, in which S(t) flips between the two states 
S = 0 and S = 1 at random times. The statistics of 
these transitions is governed by the time-dependent 
rates V(t)/r and R(t), both of which have units of 
inverse years. The underlying biophysical elements of 
AF initiation and termination we wish to describe (see 
Sec. 2.1) are encapsulated in the detailed mathemat¬ 
ical forms of these rates. It is important to keep in 
mind that these depend on the prior history of the 
patient, and so V(t) and R(t) are not set quantities, 
but random variables themselves. We will describe 
the mathematical forms of these rates next. In doing 
this we will first define the general mathematics, and 
not specify the numerical values of the various model 
parameters we are about to introduce. These will be 
discussed further below. 


3.2 Vulnerability V(t) 

The vulnerability V(t) is made up of three compo¬ 
nents; genetic predisposition, age and co-morbidity re¬ 
lated vulnerability, and episode-induced vulnerability 
[12]. We write this as 

V t = Vo + Kge + Kpi- (1) 

The quantity Vo is the vulnerability associated with ge¬ 
netic predisposition, it is time-invariant and may vary 
from patient to patient. The second component, V age , 
is the vulnerability associated with the age and with 
co-morbidities that may develop with time. In the 
present study we do not model co-morbidities inde¬ 
pendently or explicitly, and simply assume that V age 
depends primarily on age (i.e. other co-morbidities are 
implicitly correlated with age), but not on the prior 
AF-history of the subject. This type of vulnerability 
increases with time. More specifically we assume a 
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Age/co-morbidity-dependent vulnerability is low in 
the initial phases of a subject’s trajectory (i.e. t ~ 0), 
and it then increases to plateau, T4, at higher ages. 
The parameter t c describes the location of the turning 
point in this evolution, and td characterises the width 
of the sigmoid, i.e. the duration of the transition be¬ 
tween low vulnerability and the final plateau. 

It remains to define Wpi as the vulnerability associ¬ 
ated with acute remodelling. The following piecewise 
process is adopted to model V ep {. When the patient 
is not in AF, i.e. at times when S(t) = 0, we assume 
that Wpi decays exponentially with rate /?, 

(t) = -/3V epi (t), when S(t) = 0. (3) 

When a patient is in AF (i.e., when S(t) = 1), we 
assume that the vulnerability grows exponentially in 
time limited by a maximum vulnerability of V max . 
Mathematically we write 


jV ev i(t) = a [V max - Vpi(*)], When S(t) = 1. (4) 

The model parameters a and /? are the relaxation rates 
of Wpi back to 0 when S(t) = 0, or to maximum value 
when S(t) = 1. 

We further assert that the vulnerability V(t) is con¬ 
tinuous: after a random switching event into AF or 
out of AF, the initial condition of V ep { (t) in the fol¬ 
lowing segment is set as the final value of V ep i(t) right 
before the switching event. A schematic diagram of 
the evolution of V(t) is presented in Fig. 2. 


3.3 Recovery rate R(t) 

Next, we define the recovery rate R(t) to consist of 
two components: 

Rt •= Rage (t) + Repi{t). (5) 

The variable R age (t) stands for a naturally declining 
ability to recover from an established episode. The 
decline is due to age, and it is taken to be exponential 
with rate A, 

-Rage(i) := Roe~ Xt - (6) 

The second component, R e p i, is the history-dependent 
recovery rate, which is known to be a function of pre¬ 
ceding episodes of AF [11]. Similar to the above proce¬ 
dure for the vulnerability we define a piecewise process 
for R e P i. When the subject is not in AF, we assume 

= -vRepiit) when S(t) = 0, (7) 

which describes the recovery from future episodes of 
AF, and indicates an exponential decay of R epi with 
rate v during times the patient is in SR. 

When the patient is in the AF state (S(t) = 1) we 
assume that the total recovery rate, R(t) = R age (t) + 


R e p i(t) falls exponentially with rate /a 

~^R(t) = when S(t) = 1. (8) 

This different formulation is intended to capture the 
short term remodelling that takes place during an 
episode of AF, and acts to reduce the likelihood of ter¬ 
mination. Crucially we assume fi > A, i.e. the drop in 
total recovery rate during AF is quicker than the drop 
of R age (t) when the subject is in sinus rhythm. Typ¬ 
ically the total recovery rate R(t) will fall below the 
baseline rate i2 age (t), hence R ep i(t) = R(t) — R age (t ) 
formally turns negative; we comment on this below. 
The dynamics illustrated in Fig. 2. 

In order to capture the clustering of several short 
episodes frequently observed clinically after a long AF 
episode, we introduce an additional ‘boost’ term into 
the dynamics of the recovery rate. Specifically we use 
the following procedure when an episode of AF ter¬ 
minates. Say the termination time is t. Prior to this 
moment the patient will have been in AF, and so R(t) 
is governed by Eq.(8). The age-related contribution is 
Rage (t) = Roe~ xt , see Eq.(6). Together these equa¬ 
tions define R ep i = R — R age just before the end of 
the episode at t. After a long episode of AF this net 
rate R e p i will often be negative. To capture the clus¬ 
tering we assume that immediately after the end of 
an episode the total recovery rate R(t) overshoots and 
exceeds the natural background R age {t). In order to 
implement this we introduce an instantaneous boost 
at the end of an episode: we change the sign of R ep i{t) 
and we allow for an additional factor B > 0 in its mag¬ 
nitude: R e pi —> — B x R e pi- This is only implemented 
if R ep i(t) < 0 just before the end of the episode, other¬ 
wise no reflection or boost are applied. We note that 
negative rates of R e p i do not obviously have a direct 
meaning in isolation. In our model only the total re¬ 
covery rate R = R age + R e p i, has a direct biological 
interpretation. This rate is always positive. 

The different features of the dynamics are illustrated 
in Fig. 2, where we also show stylised possible time 
courses of the vulnerability and recovery rate. 


3.4 Implementation 

The model constitutes a continuous-time stochastic 
process with discrete and continuous variables. During 
episodes and during times of sinus rhythm the rates 
R(t) and V(t) evolve deterministically according the 
dynamics specified by the ordinary differential equa¬ 
tions above. The switching dynamics between 5 = 0 
and 5 = 1 (SR and AF) is stochastic and governed 
by these rates. The transition rates evolve determin¬ 
istically over time conditional on the current state of 
the subject. The model can be implemented relatively 
straightforwardly using a modified Gillespie algorithm 
[29, 30] to account for the time-dependent rates gov¬ 
erning the statistics of the switching events into and 
out of AF. This is discussed further in the Appendix. 
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Parameter 

Description 

Value 

Unit 

Comment/Reference 

1/t 

Triggering rate 

1 

1/year 

Basic simulation time unit 

Vo 

Gene, predis. vul. 

10 -t> 

n/a 

Set arbitrarily small 

Vi 

Max. age/disease vulnerability 

2920 

n/a 

Paroxysmal AF patients have up to 8 episodes/day [31] 

Vmax 

Max. episode vulnerability 

2 x vi 

n/a 

AF doubles VW in tissue [32] 

tc 

Age/disease transition time 

70 

years 

Mean incidence rate [33, Table 1] 

td 

Age/disease transition width 

3 

n/a 

Gives range ±10 years 

a 

Episode vul. saturation rate 

52 

l/yr 

AERP stabilises after one week [34] 

p 

Episode vul. relaxation rate 

182.5 

1/yr 

AERP recovers after two days [35] 

Ro 

Peak recovery rate (r.r.) (of new born) 

0.5 

second 

Min possible period of AF 

X 

Nat. degradation rate of r.r. 

1.2/50 * log(840) 


Set so that patient has 14.4 min episode at age 50 

[i 

Degradation rate of total r.r. in AF 

a 


Follows a 

V 

Relaxing rate of epi. r.r. when Healthy 

p 


Follows /3 

B 

Boost factor of r.r. post AF 

1 

n/a 

1: overshootxl; 0: no o/s 


Table 1: The parameter set used in model simulations. Where possible, parameters were estimated from literature; 
certain parameters were set to reproduce incidence behaviours from epidemiological studies. 


3.5 Model parameters 

The parameter set used in the present study is listed in 
Table 1 and discussed below. We stress that these are 
not definitive population-representative values, and 
that further exploration of the parameter space is nec¬ 
essary. 

The external triggering rate 1/r scales the vulner¬ 
ability V{t) and has units of inverse time. We here 
use r = 1 /yr, and use years as our unit of time. In 
the following we will occasionally omit units in activa¬ 
tion and recovery rates, as they are always understood 
to be per year. The numerical values of the different 
components of V{t) can then be set individually, as 
shown in Table 1. The parameters Vo, V\ and Wiax 
are dimensionless quantities. We are focused on inves¬ 
tigating effects of natural history and of AF-dependent 
remodelling, and thus set the vulnerability due to ge¬ 
netic predisposition Vo to be very small at 10 -5 . Age¬ 
ing is a key predictor of AF vulnerability [36], with in¬ 
cidence and prevalence rates increasing with age. For 
the vulnerability due to ageing/co-morbidity, V a ge (t), 
we set the maximum value V\ as 365 x8 = 2920 on the 
basis of Hoffmann et al [31] who reported a median of 
8 episodes per day in a patient cohort with paroxysmal 
AF. For vulnerability due to AF history W p i, we cite 
a recent computational study from Colman et al [32], 
who predicted that chronic AF doubles the vulnerable 
time window (in which the atria is susceptible to for¬ 
mation of AF-inducing re-entrant circuits) compared 
to baseline; we take this to assume that Vmax = 2 x V\. 

To identify the parameters of the sigmoid curve de¬ 
scribing increase in vulnerability due to ageing, we 
looked at the incidence rates of AF in several pop¬ 
ulation studies [33, 37, 38, 39] and chose parameters 
t c = 70 and G = 3, representing a sigmoid centred 
at 70 years of age, with a range of approximately 10 
years. 

The saturation and recovery rates for episode- 
related vulnerability were estimated from a series 
of experimental studies establishing the existence of 
AF-induced AF over timescales of days to weeks 
[35, 11, 24]. They analysed the atrial effective refrac¬ 
tory period (AERP), a clinical biomarker of atrial vul¬ 
nerability, over periods of artificially induced AF and 
looked at changes relative to baseline along with recov¬ 
ery rates post pacing. AERP is known to increase in 
the right atria with age [40] but shortens when in AF. 


From these studies, we set <a, the episode vulnerabil¬ 
ity saturation rate, to be 52, corresponding to a time 
period of 1 week [34] and /3, the episode vulnerability 
relaxation rate to 182.5, corresponding to a maximum 
relaxation period of 2 days [35]. 

For the parameters governing the recovery function 
R(t), we chose to estimate these based on intuition of 
the AF remodelling process. We set the peak recov¬ 
ery rate po as 0.5 seconds, which we assume to be the 
shortest possible duration of an AF episode - this cor¬ 
responds to the fastest possible recovery time from AF 
of a (new born) patient in perfect health. The decay 
rate of the natural AF recovery function, A was de¬ 
rived on the basis that, from R age = 0.5at£ = 0,a 
typical episode of a 50 year old patient would be on 
the order of 14.4 minutes in duration (1% of a day). 
This parameter would then be A = log(840 x 2)/50; 
in our simulation, we used a slightly larger value of 
A = 1.2 x log(840)/50. 

Finally, for the episode-dependent recovery rate 
R e pi (t), we assume that the degradation rate p follows 
the saturation rate of episode-dependent vulnerability, 
<a, and that the relaxation rate v similarly follows the 
corresponding relaxation, f3 in the vulnerability func¬ 
tion. For the boost factor B (the overshoot above base¬ 
line immediately after an episode terminates), we set 
this to 1, which exactly reflects the decreased ampli¬ 
tude in recovery rate from baseline due to the present 
episode of AF. 

4 Model Outputs 
4.1 General remarks 

The primary output of simulations of the model dy¬ 
namics are patient trajectories. Each run generates a 
series of AF episodes over several decades; typically 
we see gradual progression of AF from paroxysmal to 
permanent. These synthetic data differ from existing 
clinical AF data: real-world AF ECGs are sampled at 
several kilohertz over short monitoring periods from 
hours to months, whilst our data simplifies the ECG 
in to sequences of binary AF/sinus rhythm episodes. 
The theoretical resolution of the model is dictated by 
the shortest episode generated by the model parame¬ 
ters. 

We stress that the model parameters used are the 
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Figure 3: Time series of the progression. A sample path is visualized by the red line, and from top to bottom: 
(a) Binary signal whether the patient is in an AF episode, (b) the fraction of the time in AF episodes per day, 
(c) the fraction of the time in AF episode per week, (d) the fraction of the time in AF episode per month. In 
(b-d) 9 other sample paths are plotted in grey lines in order to illustrate the distribution induced by intrinsic 
stochasticity. Inset: the progression of AF episode frequency and duration at time points in the simulation. 


same for all simulated patients. Each virtual patient 
has characteristics identical to the population average. 
Thus, our model does not capture heterogeneity in 
those parameters. A systematic analysis of the full pa¬ 
rameter space is beyond the scope of the present work, 
but below we highlight a number of data which hint 
at the range of possible trajectories that this model 
generated. It should be noted that while the parame¬ 
ters as listed in Table 1 are the same for all subjects 
in the model, different runs of the simulation will lead 
to different sequences of AF/SR episodes. This het¬ 
erogeneity comes from the inherent randomness of the 
initiation and termination dynamics in the model. 

Results reported below are from simulations of 5000 
independent patients. On each simulation run, the 
patient lives to age 100 years, develops AF at some 
point and is found to enter permanent AF in their 
lifetime. This is not an accurate reflection of reality, 
of course, but this allows us to analyse the different 
patient trajectories as a consequence of both long- and 
short-term remodelling. 

4.2 Patient trajectories and AF pro¬ 
gression 

The top left panel of Fig. 3 illustrates a sample time 
series plotted over a 30-year time span. For any point 
in time the patient is either in AF or in sinus rhythm; 
episodes of AF are highlighted in red. In this example, 
an initial period of quiescence gives way to a sequence 
of AF episodes of increasing frequency and duration. 
As can be seen on the magnification on the right-hand 
side of the figure, where each time series represents 
a single month, there are many short and clustered 
episodes which may not be visible over the longer time- 
scale view. These eventually become longer in dura¬ 
tion until finally the patient remains in AF constantly 
until death (not shown). 

The red (dark) curves in panels (b)-(d) of Fig. 3 we 
show the daily, weekly, and monthly AF burden of the 
patient trajectory shown in panel (a). Following stan¬ 
dard conventions [31, 41], burden is here defined as 


the fraction of time spent in AF during a given obser¬ 
vation window. To illustrate the variability of model 
outcomes panels (b)-(d) depict 9 other sample paths 
(light grey). AF burden shows a monotonic increasing 
trend. Naturally the patient-to-patient variance re¬ 
duces for longer observation windows (monthly versus 
daily). 

Using the parameter set in Table 1, we find that all 
patients develop paroxysmal AF at some point, and 
that they eventually enter permanent AF before the 
endpoint of the simulation at age 100. 

The model then allows us to simulate population 
statistics for the age at which patients develop parox¬ 
ysmal AF (‘event 1’) and for the time at which they 
transition into permanent AF (‘event 2’). These are 
shown in Fig. 4 (panels (a) and (b)), along with 
statistics of the time that elapses between the on¬ 
set of paroxysmal and permanent AF (panel (c)). It 
seems reasonable to assign a minimum threshold which 
would prompt a patient to present themselves for clin¬ 
ical diagnosis (assuming symptomatic AF). We set 
this threshold as an AF episode of duration greater 
than 14.4 minutes (1% of a day). The model results 
show that the first episode of AF lasting longer than 
the threshold typically occurs when patients are aged 
around 44-48 years, which may be perceived as earlier 
than recorded in the literature [33, 37, 38, 39]. Patients 
subsequently progress to permanent AF over a period 
of 28-32 years (panel (c)). A potential interpretation 
of these results is discussed in Section 5. 

4.3 Behaviour after onset of paroxys¬ 
mal AF 

We now turn to a discussion of model outputs in a 
form that may be clinically relevant. In our model we 
assume that the starting point of a patient’s clinical 
trajectory is the time at which they develop parox¬ 
ysmal AF. Given the intrinsic stochasticity, this time 
point will vary from patient to patient, as indicated in 
panel (a) of Fig. 4. 

For any one patient we can use the time at which 
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Figure 4: Distributions of age at which paroxysmal AF 
sets it (panel (a)), age at which permanent AF sets in 
(panel (b)), and the time elapsed between these two 
events (panel (c)). Data were generated from 5000 
independent simulation runs of the model. 


they enter paroxysmal AF as a reference time to follow 
their future progression. Data are shown in 5; the 
horizontal axes in the top panels of the figure indicate 
time after the onset of paroxysmal AF. 

In the top left panel of Fig. 5, annual AF burden 
is shown. The burden increases from near zero to 
one, representing the increase in the proportion of time 
spent in AF as the condition progresses. Each point 
on the curve represents the mean burden over 5000 
simulated patients. 

In the lower-left part of the figure we shown his¬ 
tograms of AF burden at the 15 time points indicated 
in the top-left panel. The numbers in the smaller 
panels refer to the number of years since the onset 
of paroxysmal AF. At initial times after the onset of 
paroxysmal AF the burden shows a unimodal distri¬ 
bution, and over the following 17 years the burden in¬ 
creases for some patients while remaining low for oth¬ 
ers producing a left skewed distribution. Around year 
19 there is a shift towards a higher burden and the dis¬ 
tribution becomes more uniform, indicating that those 
patients that previously had low burden are now ex¬ 
periencing the effects of remodelling. By year 21these 
effects are evident in even more patients and the dis¬ 
tribution becomes right skewed. Finally over years 
24 — 37 all patients experience the effects of remod¬ 
elling and progress to permanent AF indicated by the 
distribution tending to a burden of 1. 

The right-hand side of Fig. 5 shows the total num¬ 
ber of episodes per year (top panel) for the same 
dataset, along with detailed breakdowns of the pop¬ 
ulation statistics in the lower right-hand panels. The 
number of episodes per year rises sharply and peaks at 
around 20 years after the onset of paroxysmal AF, and 
then gradually decreases. This reflects the observation 
that AF episodes begin intermittently, and become 
more frequent as structural remodelling increases vul¬ 
nerability and decreases recovery rate. Consequently 
the episode durations lengthen and the inter-episode 
times shorten, so the number of episodes is reduced 
though the time spent in AF (AF burden) increases 
(left-hand panels of Fig. 5). 

Initially the number of episodes experienced per year 
follows a unimodal distribution, with a peak at low 
numbers. The pace with which AF progresses in differ¬ 
ent virtual patients varies considerably; some will still 
have infrequent episodes so a low number of episodes 


per year, others will have experienced a moderate 
degree of re-modelling and will have many very fre¬ 
quent episodes. A third group will have experienced 
substantial structural re-modelling and so they have 
fewer episodes of longer duration. This leads to the 
broad distribution at around 17 years after the onset of 
paroxysmal AF. Eventually all patients will progress 
to permanent AF and therefore have fewer episodes 
(19-27 years after threshold). Finally AF is perma¬ 
nent, i.e. there is one continuous episode, and then 
the centre of the distribution moves towards 1 single 
episode per year. 

4.4 Behaviour prior to onset of perma¬ 
nent AF 

We now change perspective and look at patient pro¬ 
gression relative to the point in time at which they de¬ 
velop permanent AF. As discussed, different patients 
can progress differently after the onset of paroxysmal 
AF, and there will be considerable variation of the time 
at which they enter permanent AF. To allow for the 
identification of potential patterns that can determine 
progression to permanent AF, we have replotted the 
model outputs in Fig. 6. Time is now defined relative 
to the point at which permanent AF sets in, the hori¬ 
zontal axis in the top panels now depict time until the 
onset of permanent AF. The left-hand panel shows the 
annual AF burden, for the time leading up to the mo¬ 
ment of transition to permanent AF, the panels on the 
right show data for the number of episodes experienced 
per year. 

From the left-hand panels in Fig. 6 we can see that 
the annual burden begins to increase very slowly until 
16 years prior to transition to permanent AF. At that 
point, patients start to progress at different speeds 
causing the distribution to become more unimodal 
(years 14 — 12 prior to transition). By 10 years before 
transition the patients are experiencing fairly high bur¬ 
den, which continues to gradually increase until they 
enter permanent AF. 

Examining the right-hand panels in 6 one sees that 
the number of episodes at 17 years prior to transition 
has noticeably increased but the corresponding burden 
(left-hand panel) is still low. This indicates that there 
are many short episodes - and we know that these are 
particularly unlikely to be picked up clinically simply 
due to timing. The number of episodes continues to 
increase and at 14 years prior to transition the dis¬ 
tribution for the number of episodes becomes more 
unimodal and symmetrical (ignoring a peak at zero 
episodes, reflecting patients who have not yet entered 
AF). The distribution shows a wide range in the num¬ 
ber of episodes, with the corresponding burden also 
starting to increase. This illustrates that some people 
are experiencing the effects of short-term remodelling, 
and having many clustered episodes. At 12 years be¬ 
fore transition we see that most patients have many 
episodes, but burden still follows a relatively normal 
distribution. This indicates that some patients are ex¬ 
periencing the effects of long-term re-modelling and 
the duration of the individual episodes are increasing, 
while other patients are only experiencing the effects 
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Figure 5: Left: Top panel shows annual AF burden after patient has entered paroxysmal AF. Horizontal axis 
indicates time elapsed since onset of paroxysmal AF. The mean of 5000 simulation runs is indicated as solid line, 
grey area shows the standard deviation. Smaller panels at the lower left show distribution of burdens across the 
patient population at the indicated points in time (years after onset of paroxysmal AF). Right: Analogous data 
for the number of AF episodes experienced in a moving window of 12-months duration. 


of short-term re-modelling and although they experi¬ 
ence many episodes they are of short duration. From 
6 years until the transition to permanent the number 
of episodes decreases as the duration increases, until 
eventually the burden is 1 and the patient remains in 
AF. 

5 Discussion 

In this study we have presented an individual-based 
model of AF progression. The model operates at pa¬ 
tient level and operates a coarse-grained description of 
their progression through the different stages of AF. 
The fundamental variable of the model is the current 
AF status of the patient, being either in AF or in si¬ 
nus rhythm. The transitions between these states are 
governed by history-dependent activation and termi¬ 
nation rates, respectively. These rates are constructed 
to retain some key features of the physiology and bio¬ 
physics of AF; in particular activation and termina¬ 
tion rates both include terms describing acute elec¬ 
trical and structural remodelling. If a patient has 
a history of AF they are more likely to initiate fur¬ 
ther episodes, and these have a reduced probability 
to terminate. Other factors we include are a baseline 
vulnerability and age-dependent factors. The compo¬ 
nents related to age include co-morbidities which may 
develop over time, although we do not model these 
in detail. Our model builds on the conceptual frame¬ 
work proposed by Heijman et al [12] and we remark 


that this conceptual model discusses only AF initia¬ 
tion and maintenance, but does not explicitly detail 
AF termination. 

The formulation of our model is based on reaction 
kinetics with transitions between states that are mod¬ 
ulated by the activation and recovery functions. This 
type of approach has been widely used in other areas 
(for example to describe the state of ion channels in 
cardiac cells), and yields behaviours that are broadly 
consistent with experimental observation. In the case 
of our model, we argue that it is reasonable to as¬ 
sume that patients with AF will transition between 
SR and AF states, and that the transition rates are 
affected by age and remodelling. The functional form 
of our age/co-morbidity-related vulnerability is con¬ 
sistent with the incidence curves reported in clinical 
studies [33, 38]. Where possible, we have chosen pa¬ 
rameters based on evidence from the literature (Table 
1). However, in a real patient we would expect some 
parameters such as the triggering rate to be variable 
rather than fixed as in the current model. We have also 
implicitly included the effects of co-morbidities within 
our description of vulnerability, which may not nec¬ 
essarily be age-correlated. Separating out the effects 
of co-morbidities in the vulnerability would allow us 
to investigate the progression of lone AF compared to 
other patient groups [42]. We expect that future stud¬ 
ies will address these extensions of the model. 

From the model outputs, we have analysed the typ¬ 
ical disease progression of a virtual patient generated, 
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Figure 6: (Left panel) Burden before patient goes into the permanent episode. (Right panel) Number of episodes 
before the patient goes into the permanent episode. In both panels we present the mean of the random variables 
(solid black) and mean db 1 standard deviation (grey area) at the top, and at the bottom we measure the 
distributions of the random variables at 15 sampled time t = [21,19,17,16,15,14,13,12,10, 8, 7, 6, 5, 3,1] years 
before the permanent episode. 


and we have characterised the statistics derived from 
an idealised cohort of this population-averaged pa¬ 
tient. In particular we report data relating to AF bur¬ 
den and the number of episodes experienced per year 
from the onset of paroxysmal AF forward. Secondly we 
have studied population statistics during the build-up 
to developing permanent AF. We did not study persis¬ 
tent AF in detail within the present work, other than 
to define the diagnosis threshold; the ability to predict 
progression to persistent AF would be of further clin¬ 
ical interest, given the difficulty of making accurate 
diagnoses from limited clinical monitoring windows. 

The results from the model appear to show that 
the first episode with a duration greater than 14.4 
minutes, our diagnosis threshold for onset of parox¬ 
ysmal AF, starts earlier in the patient’s life than re¬ 
ported in the literature [33, 38]. However, it must be 
noted that the model records all events, while in real¬ 
ity these events may be asymptomatic, and therefore 
remain undiagnosed and therefore undocumented for 
many years [43]. As a consequence of the recording 
of asymptomatic events prior to clinical diagnosis, the 
time taken to transition to permanent AF (Fig. 4(c)) 
may appear to be overly long. The follow-up period 
in published studies on the transition to permanent 
AF (e.g. [44, 45]) are of a duration of between 1 and 5 
years or up to 30 years in studies of lone-AF (e.g. [42]), 
and include patients in both paroxysmal and persistent 
states at recruitment, so direct comparison with the 
present model result is challenging. However if we con¬ 


sider the age of transition of the patient to permanent 
AF, see Fig. 4(b), this does reflect findings in the liter¬ 
ature that age is an independent predictor of chronic 
AF (both persistent and permanent) and that the risk 
associated with age increases beyond the age of 74 [5]. 
Additionally, this model takes no account of any med¬ 
ical interventions where the transition to permanent 
AF may be affected by different treatment strategies. 
Combining these results with models, such as Lord et 
al. [10], that capture the longer-term impacts of med¬ 
ical interventions on a healthcare system may provide 
insights into when different treatment strategies are 
best applied along the progression trajectory. 

We stress that the model as it stands is of a stylised 
nature. While model parameters are chosen in agree¬ 
ment with existing clinical literature, we have made 
no attempt at validating the model outcome against 
real-world patient data. We do believe however that 
models of the type we propose here are placed some¬ 
where between detailed mechanistic approaches and 
large-scale epidemiological studies, provide a promis¬ 
ing way forward in integrating modelling of AF. We 
also note that analogous approaches are indeed possi¬ 
ble for other diseases which progress through different 
stages [46, 47, 48]. 

Now that a basic model is in place future work may 
focus on extending and refining it. Elements to in¬ 
corporate in the future may include a more detailed 
breakdown of AF triggers, patient vulnerability and 
the recovery rate. This would ultimately have to be in- 
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formed by the outputs of mechanistic models. Further 
important aspects include patient-to-patient variation 
of the underlying parameters, and ultimately the mod¬ 
elling of intervention techniques such as rhythm or rate 
control [49] At the other end, comparison and valida¬ 
tion against long-term clinical cohort data would be 
desirable. 

If progress can be made along these lines, one 
may hope that the intermediate modelling approach 
we propose here will ultimately establish a link be¬ 
tween small-scale mechanistic models and large-scale 
epidemiological data. Properly validated individual- 
based models would ultimately contribute to the de¬ 
velopment of computer-based personalised medicine, 
and targeted pharmacological or surgical intervention. 
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Drawing At from this distribution is non-trivial, and 
so we adopt a trick. While the onset of an episode 
in the clinic can obviously not be attributed to one 
unique the component to P (i.e. either to Vo, Page 
or P epi ) we formally treat these as separate events in 
the simulation. I.e. we draw on waiting time, Ati to 
an episode triggered with constant rate Vo/r, another 
waiting time At 2 generated using P age /r, and another 
(Ats) based on V ep i/r. How this is done in practice 
is explained below. We then use the shortest of these 
three times, At = min{Ati, At 2 , Ats}. This faithfully 
reproduces the statistics governed by the total rate, 
i.e. that of Eq. (10). 

Simulation time is then advanced to t + At, the re¬ 
covery rates are forward integrated, and an episode is 
activated at t + At (i.e. we set 5 = 1). 

It remains to explain how to generate waiting times 
until events triggered by the above three rates. 

Given that Vo/r is constant the distribution of wait¬ 
ing times for events follows an exponential distribution 
with parameter Vo/r. This can be sampled directly 
setting Ati = — ^ In Ui, where U\ is a random vari¬ 
able drawn uniformly from the interval ( 0 , 1 ]. 

To construct At 2 , one plugs in V age (t') = 
V\/[r (1 + exp [—{t' — t c )/td\)\- Integration gives 
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A sample of At 2 can then be obtained by drawing a 
uniform random number U 2 from the unit interval, and 
by solving 


A Further details of the simulation al¬ 
gorithm 

Simulation of the model generates sample paths, char¬ 
acterised by an alternating series of episode initiation 
and termination times. We here briefly describe how 
these are obtained. 

A.l Activation of episodes 

Assume that the individual is in sinus rhythm at time 
£, i.e. S(t) = 0 . Its current vulnerability and its re¬ 
covery rate are known from the subject’s prior history. 
The next step of the simulation is now to find the time 
At until the onset of the next episode. The statistics 
of At are governed by the time evolution of P from 
time t forward. 

The subject is not in AF at time £, and so until the 
next episode sets in the time-dependence of vulnera¬ 
bility is given by (t' > t) 


P(t / ) = P0 + Page(O + Pepi(O- (9) 


The waiting time until the next episode, At, is then 
governed by the following cumulative probability dis¬ 
tribution 


Pr(At < T) = 1 — exp 



( 10 ) 


exp 


ViU 


log 


3 (t+At 2 )/£d _|_ e t c /td 

(?t/td -|- g^cAd 


= U 2 . (12) 


This can be done in explicit form. 

Finally we construct At 3 , using the same idea. We 
have Pepi(t') = [P(t) - Pi (t) - Vo] exp (-/?(£' - t)) us¬ 
ing Eq. (3), and keeping in mind that the subject is 
not in AF. After integration we see that At 3 can be 
obtained from setting 


exp 


V(t) - Pi ft) - Pp ^_p t _ p -{3(t+At 3 ) 


r/3 


’) 


= u 3 . 

(13) 


Again, this can be inverted in explicit form when 


U3 > exp 


v(t) - u(t) - Vo __ et 


(14) 


and we define At 3 00 otherwise. 


A.2 Termination of episodes 

Now assume the subject is in AF at time t. The al¬ 
gorithm proceeds by drawing a time At to recovery. 
The recovery rate for t' > t follows the exponential 
law R(t') = Rft)e _/X ^ during an episode, and so 
we can determine At < 00 from 


exp 



^g-/x£ _ e -/x(t+At)^ 


= U, 


(15) 
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when 


U > exp 


T 


(16) 


and At —>• oo otherwise, with a uniform random num¬ 
ber U drawn from the unit interval. 

Once At is drawn, activation and recovery rates are 
integrated forward until t -h At, and the episode is ter¬ 
minated (S(t +At) = 0). The algorithm then proceeds 
as described in the previous section. 
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